
-- {-# LANGUAGE ViewPatterns        #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE RebindableSyntax    #-}
-- {-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
module Main where

import Codec.Picture
import Codec.Picture.Types
import Codec.Picture.Jpg
import qualified Codec.Picture.Metadata as Meta
import Codec.Picture.Extra
import Data.Array.Accelerate.Data.Complex as A
import Data.Array.Accelerate as A 
import qualified Data.Vector as V

import Data.Time.Clock

import qualified Data.ByteString.Lazy as BL
-- import qualified Data.ByteString.Lazy.Char8 as BLC
import qualified Data.ByteString as B

import Data.Array.Accelerate.IO.Codec.BMP             
import Data.Array.Accelerate.IO.Codec.Picture
import Data.Array.Accelerate.LLVM.Native as CPU
-- import Data.Array.Accelerate.LLVM.PTX as PTX
import Data.Array.Accelerate.Data.Complex 
import System.IO

import Foreign
import Foreign.C
import Prelude
import qualified MyLib (someFunc)

import qualified Data.Array.Accelerate.Control.Lens as ACL


scale_image :: Acc (Array DIM2 PixelRGB8)-> Exp Int ->Exp Int->Acc (Array DIM2 PixelRGB8)
scale_image image_data nwidth nheight = A.generate (A.I2 nheight nwidth)
    (\ix-> inter_ (compute_conti ix ) )
    where 
      A.I2 height width = A.shape image_data 
      compute_conti::Exp DIM2 ->  (Exp Float ,Exp Float)
      compute_conti index= ( A.fromIntegral x /A.fromIntegral nwidth,A.fromIntegral y /A.fromIntegral nheight )
        where 
          (Z :. y :. x)=unlift index

      mult_pix::Exp Float-> Exp PixelRGB8 ->Exp PixelRGB8
      mult_pix m (PixelRGB8_ a1 a2 a3)  = PixelRGB8_ (A.floor ( A.fromIntegral a1 *  m )) (A.floor ( A.fromIntegral a2 *  m )) (A.floor ( A.fromIntegral a3 * m ))
    
      pix2float::Exp PixelRGB8-> (Exp Float,Exp Float,Exp Float)
      pix2float (PixelRGB8_ a1 a2 a3)= (A.fromIntegral a1,A.fromIntegral a2,A.fromIntegral a3)


      float2pix:: (Exp Float,Exp Float,Exp Float)->Exp PixelRGB8
      float2pix (a1,a2,a3)= PixelRGB8_(A.floor a1) (A.floor a2) (A.floor a3) 

      multfloat:: Exp Float ->(Exp Float,Exp Float,Exp Float)->(Exp Float,Exp Float,Exp Float)
      multfloat m (a1,a2,a3)= (m*a1,m* a2,m*a3)

      plusfloat:: (Exp Float,Exp Float,Exp Float) ->(Exp Float,Exp Float,Exp Float)->(Exp Float,Exp Float,Exp Float)
      plusfloat (b1,b2,b3) (a1,a2,a3)= (b1+a1,b2+ a2,b3+a3)

      inter_:: (Exp Float,Exp Float)-> Exp PixelRGB8
      inter_   (x,y) =
        float2pix $ plusfloat  (multfloat (lambdax*lambday) dp1) $ plusfloat (multfloat (lambdax*lambday_)  dp2) $ plusfloat   (multfloat (lambdax_*lambday ) dp3 )  (multfloat (lambdax_*lambday_ ) dp4)               
        where 
          x1=A.min  (width-1) (A.floor (A.fromIntegral width *x))
          x2=A.min  (width-1) (A.floor (A.fromIntegral width *x) + 1  )
          y1= A.min (height-1) (A.floor (A.fromIntegral height *y))
          y2= A.min (height-1) (A.floor (A.fromIntegral height *y)+1)
          dx1=A.fromIntegral x1 /A.fromIntegral width ::Exp Float
          dx2=A.fromIntegral x2 /A.fromIntegral width ::Exp Float
          dy1=A.fromIntegral y1 /A.fromIntegral height :: Exp Float
          dy2=A.fromIntegral y2 /A.fromIntegral height :: Exp Float

          lambdax=A.min 1.0  (( dx2-x)*A.fromIntegral width) 
          --lambdax=0.5
          lambdax_ = 1 - lambdax
          lambday= A.min 1.0 ((dy2-y)*A.fromIntegral height)
          --lambday=0.5;
          lambday_ = 1 - lambday
          dp1=pix2float $ image_data A.! (lift (Z :. y1 :. x1))
          dp2=pix2float $ image_data A.! (lift (Z :. y2 :. x1))
          dp3=pix2float $ image_data A.! (lift (Z :. y1 :. x2))
          dp4=pix2float $ image_data A.! (lift (Z :. y2 :. x2))




downsample_image:: Acc (Array DIM2 Float)-> Acc (Array DIM2 Float)
downsample_image  image =let 
    I2 h w = shape image 
  in
    backpermute (I2 (A.quot h 2 +1 ) (A.quot w 2 + 1)) (\(I2 x y)-> (I2 (A.min (h-1) (2*x) ) (A.min (w-1) (2*y) ) ) ) image   





gauss_convolution1 :: Acc (Array DIM1 Float) -> Acc (Array DIM2 Float) -> Acc (Array DIM2 Float)
gauss_convolution1 gauss_op image =temp_image2
  where
    A.I1 n  = A.shape gauss_op
    -- Size of the full image
    A.I2 h w = A.shape image

    myfun :: Exp DIM3 -> Exp Float
    myfun (A.I3 y x kery) =
      (gauss_op A.! (A.I1 kery)) * (image A.! (A.I2 (A.min (h - 1) (A.max 0 (y + kery - (n `A.quot` 2) ))) x))
    temp_image1= A.fold (A.+) 0 $ A.generate (A.I3 h w n ) myfun

    temp_image2= A.fold (A.+) 0 $ A.generate (A.I3 h w n ) myfun1

    myfun1 :: Exp DIM3 -> Exp Float
    myfun1 (A.I3 y x kery) =
      (gauss_op A.! (A.I1 kery)) * (temp_image1 A.!   (A.I2 y (A.min (w - 1) (A.max 0 (x + kery - (n `A.quot` 2)  ))) ) )


normalize_gauss_plat1::Acc (Array DIM1 Float)->Acc (Array DIM1 Float)
normalize_gauss_plat1 gauss_op = A.map  myfun gauss_op
  where 
    ave= A.the $ A.foldAll (A.+) 0 gauss_op
    myfun:: Exp Float ->Exp Float
    myfun x= x/ave


generate_gauss_op1:: Exp Float-> Acc (Array DIM1 Float)
generate_gauss_op1  sig= generate (lift (Z :. m )) compute_conti
  where 
    m= 2*A.ceiling (  3*sig) +1  :: Exp Int
    m1= A.fromIntegral m ::A.Exp Float
    compute_conti::Exp DIM1 -> Exp Float
    compute_conti index=  (A.exp  (-((y1-m1/2.0)**2 ) /(2*sig*sig))) / ((2*A.pi )**0.5*sig) 
      where  
        (Z :. y )=unlift index 
        y1= A.fromIntegral y +0.5 :: Exp Float


generate_gauss_op:: Exp Float-> Acc (Array DIM2 Float)
generate_gauss_op  sig= generate (lift (Z :. n :. m )) compute_conti
  where 
    m= 2*A.ceiling (  3*sig) +1  :: Exp Int
    n=m
    m1= A.fromIntegral m ::A.Exp Float
    n1= A.fromIntegral n ::A.Exp Float
    compute_conti::Exp DIM2 -> Exp Float
    compute_conti index=  (A.exp  (-((x1-m1/2.0)**2 +(y1-n1/2.0)**2 ) /(2*sig*sig))) /(2*A.pi*sig*sig) 
      where  
        (Z :. y :. x)=unlift index 
        x1=A.fromIntegral x +0.5 :: Exp Float
        y1= A.fromIntegral y +0.5 :: Exp Float


compute_base_gray_image::Acc(Array DIM2 Float)->Acc (Array DIM2 Float)
compute_base_gray_image image=gauss_convolution1 gauss_op image
  where 
    sig0=1.2489996
    gauss_op=normalize_gauss_plat1 $ generate_gauss_op1 sig0
    -- gauss_op=normalize_gauss_plat $ generate_gauss_op sig0


-- (A.replicate (A.lift (Z :. All :. All :. (1::Exp Int) )) image)
generate_one_octave:: Acc (Array DIM2 Float)->Acc (Array DIM3 Float)
generate_one_octave image =  myfun1 (myfun 5)  $myfun1 (myfun 4)  $  myfun1 (myfun 3)   $  myfun1 (myfun 2)  $  
  myfun1 (myfun 1 )  (A.replicate (A.lift (Z :. All :. All :. (1::Exp Int) )) image)     
  where 
    sig0=1.6::Exp Float--- 1.2489996
    k= (A.constant 2.0) **( 1/3.0)::Exp Float
    myfun :: Exp Float-> Acc(Array DIM2 Float)
    myfun i = gauss_convolution1 gauss_op image 
      where 
        gauss_op=normalize_gauss_plat1 $ generate_gauss_op1 ( ((k**(2*i))-1 )**0.5 *sig0  )

    myfun1:: Acc (Array DIM2 Float)->Acc (Array DIM3 Float)->Acc (Array DIM3 Float)
    myfun1  a2 a1=  a1 A.++ (A.replicate (A.lift (Z :. All :. All :. (1::Exp Int)) ) a2)



myifThenElse::Bool ->a->a->a
myifThenElse True x _ =x
myifThenElse False _ y =y


generate_gauss_pyramid:: Array DIM2 Float->Int->[Array DIM3 Float]
generate_gauss_pyramid image num_oc=myfun  [CPU.run $generate_one_octave$use image] 1
  where
    myfun :: [Array DIM3 Float]->Int->[Array DIM3 Float]
    myfun (x:xs) num1= myifThenElse (num1 Prelude.>= num_oc) (x:xs) (myfun ( myfun1 x:(x:xs))  (num1+1) )

    myfun1::Array DIM3 Float->Array DIM3 Float
    myfun1  octave1 = CPU.run$ generate_one_octave  ( downsample_image $ A.slice (use octave1) (A.lift (Z :. All :. All :. (3::Exp Int))) )  


generate_dog_pyramid:: [Array DIM3 Float]->[Array DIM3 Float]
generate_dog_pyramid octaves =Prelude.map myfun octaves
  where 
    myfun :: Array DIM3 Float->Array DIM3 Float
    myfun octave =CPU.run $ A.zipWith (-)  (A.init$ use octave )   (A.tail $ use octave)  



get_maxvalue::Acc (Array DIM3 Float)->Exp Int ->Exp Int ->Exp Int->Exp Float
get_maxvalue image x y z = re
  where 
    mytool::Exp Int->Exp Int ->Exp Int->Exp Float
    mytool x1 y1 z1= image A.! (A.I3 x1 y1 z1)

    mytool12:: Exp Int ->Exp Int ->Exp Float
    mytool12 x1 y1 = A.max (mytool x1 y1 (z-1) ) $ A.max (mytool x1 y1 z) (mytool x1 y1 (z+1) )

    mytool1::Exp Int->Exp Float
    mytool1 x1= A.max (mytool12 x1 (y-1)) $ A.max (mytool12 x1 y) (mytool12 x1 (y+1))

    re=A.max (mytool1 (x-1)) $A.max (mytool1 x) (mytool1 (x+1))


get_minvalue::Acc (Array DIM3 Float)->Exp Int ->Exp Int ->Exp Int->Exp Float
get_minvalue image x y z = re
  where 
    mytool::Exp Int->Exp Int ->Exp Int->Exp Float
    mytool x1 y1 z1= image A.! (A.I3 x1 y1 z1)

    mytool12:: Exp Int ->Exp Int ->Exp Float
    mytool12 x1 y1 = A.min (mytool x1 y1 (z-1) ) $ A.min (mytool x1 y1 z) (mytool x1 y1 (z+1) )

    mytool1::Exp Int->Exp Float
    mytool1 x1= A.min (mytool12 x1 (y-1)) $ A.min (mytool12 x1 y) (mytool12 x1 (y+1))

    re=A.min (mytool1 (x-1)) $A.min (mytool1 x) (mytool1 (x+1))


compute_extremum_points::Acc (Array DIM3 Float)->Acc (Vector (Int,Int,Int))
compute_extremum_points image=A.afst $ A.filter myfun1 indexs
  where
    
    A.I3 h w len=  A.shape image 
    indexs= A.generate (A.I3 (h-2) (w-2) (len-2) ) (\(A.I3 i1 i2 i3)-> lift ( i1+1 ,i2+1 ,(i3+1) ))  
  
    maxdata= A.map (\(A.T3 x y z)-> get_maxvalue image x y z ) indexs
    mindata= A.map (\(A.T3 x y z)-> get_minvalue image x y z ) indexs


    threshold=0.5*0.04/3.0 ::Exp Float

    myfun1::Exp (Int ,Int ,Int)->Exp Bool
    myfun1  x =  if (A.abs (image A.! (A.I3 i1 i2 i3 ))) A.<= threshold
      then 
        A.constant A.False
      else if image A.! (A.I3 i1 i2 i3 ) A.== maxdata A.! (A.I3 (i1-1) (i2-1) (i3-1))
      then A.constant A.True  
      else if image A.! (A.I3 i1 i2 i3 ) A.== mindata A.! (A.I3 (i1-1) (i2-1) (i3-1)) 
      then A.constant A.True 
      else 
        A.constant A.False
      where 
        (i1,i2,i3)=unlift x :: (Exp Int, Exp Int,Exp Int) 







compute_accurate_extrum_pos:: Acc (Array DIM3 Float)->Acc (Vector (Int,Int,Int))-> Acc (Vector (Float,Float,Float))
compute_accurate_extrum_pos dog extrum = A.afst $ A.filter mytool_panduan $A.map final_fun $ A.map (myf2i.myfun1.myf2i.myfun1.myf2i.myfun1.myf2i.myfun1) extrum
  where
    mytool_panduan :: Exp (Float,Float ,Float)->Exp Bool
    mytool_panduan (A.T3 xs1 xs2 xs3) = xs1 A.>0.5 A.&& xs2 A.> 0.5 A.&& xs3 A.> 0.5
     
    
    A.I3 h w len=  A.shape dog 

    myf2i :: Exp (Float,Float ,Float) ->Exp (Int,Int,Int)
    myf2i xs= lift $ (A.round i1 ,A.round i2,A.round i3) 
      where 
        (i1 ,i2,i3)= unlift xs ::(Exp Float,Exp Float,Exp Float)

    myi2f:: Exp (Int,Int,Int)->Exp (Float,Float,Float)
    myi2f xs=lift $ (A.fromIntegral i1 ,A.fromIntegral i2,A.fromIntegral i3)  
      where 
        (i1 ,i2,i3)= unlift xs ::(Exp Int,Exp Int,Exp Int)

    mytool:: Exp Int ->Exp Int ->Exp Int->Exp Float
    mytool i1 i2 i3= dog A.! (A.I3 i1 i2 i3)

    myget_d:: Exp (Int,Int,Int )-> (Exp Float,Exp Float,Exp Float)
    myget_d xs =  ( (mytool (i1+1) i2 i3 - mytool (i1-1) i2 i3 )/2.0 , (mytool i1 (i2+1) i3 - mytool i1 (i2-1) i3 )/2.0 , (mytool i1 i2 (i3+1) -  mytool i1 i2 (i3-1)  )/2.0 )    
      where 
        (i1,i2,i3)=unlift xs ::(Exp Int,Exp Int,Exp Int)
-- [Dxx,Dxy, Dxz,Dyy,Dyz,Dzz]
    myget_dd:: Exp (Int,Int,Int)->Exp (Float ,Float,Float,Float,Float,Float)
    myget_dd xs = lift $ (mytool (i1+1) i2 i3 + mytool (i1-1) i2 i3 -2 * mytool i1 i2 i3, 
      (mytool (i1+1) (i2+1) i3 + mytool (i1-1) (i2-1) i3 - mytool (i1+1) (i2-1) i3 - mytool (i1-1) (i2+1) i3)/4.0 ,
      (mytool (i1+1) i2 (i3+1) +mytool (i1-1) i2 (i3-1) - mytool (i1-1) i2 (i3+1) -mytool (i1+1) i2 (i3-1) )/4.0 ,
      mytool i1 (i2+1) i3 + mytool i1 (i2-1) i3 -2 * mytool i1 i2 i3,
      (mytool i1 (i2+1) (i3+1) + mytool i1 (i2-1) (i3-1) - mytool i1 (i2-1) (i3+1) - mytool i1 (i2+1) (i3-1) )/4.0,
      mytool i1 i2 (i3+1) + mytool i1 i2 (i3-1) -2 * mytool i1 i2 i3 )
      where 
        (i1,i2,i3)=unlift xs ::(Exp Int,Exp Int,Exp Int) 

    myget_det:: Exp (Float ,Float ,Float ,Float,Float,Float)->Exp Float
    myget_det xs=dd1*dd4*dd6+ 2*dd2*dd3*dd5-dd3*dd3*dd4-dd1*dd5*dd5-dd6*dd2*dd2  
      where 
        (dd1,dd2,dd3,dd4,dd5,dd6)=unlift xs:: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float,Exp Float)

    myget_inverse:: Exp (Float ,Float ,Float ,Float,Float,Float)-> (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float,Exp Float)
    myget_inverse xs=((dd4*dd6-dd5*dd5)/det, (dd3*dd5-dd2*dd6)/det,(dd2*dd5-dd4*dd3)/det,
      (dd1*dd6-dd3*dd3)/det,
      (dd2*dd3-dd1*dd5)/det,
      (dd1*dd4-dd2*dd2)/det ) 
      where 
        (dd1,dd2,dd3,dd4,dd5,dd6)=unlift xs:: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float,Exp Float)
        det=myget_det xs

    myget_cha::Exp (Int,Int,Int)->(Exp Float ,Exp Float ,Exp Float)
    myget_cha xs = (-d1* in1-d2*in2-d3*in3,
      -d1*in2-d2*in4-d3*in5,
      -d1*in3-d2*in5-d3*in6 )
      where 
        (d1,d2,d3)= myget_d xs 
        dd=myget_dd xs
        (in1,in2,in3,in4,in5,in6)= myget_inverse dd

    image_border_width=3::Exp Int

    -- 判断边界
    mypanduan::Exp (Int,Int,Int)->Exp Bool
    mypanduan xs= if i1 A.< image_border_width A.|| i1 A.>= (h-image_border_width) A.|| i2 A.< image_border_width A.|| i2 A.>= (w-image_border_width) A.|| i3 A.<1 A.|| i3 A.>= (len-1)
      then A.constant A.False
      else A.constant A.True
      where
        (i1,i2,i3)=unlift xs ::(Exp Int,Exp Int,Exp Int)
    
    mybound:: Exp Float->Exp Float
    mybound x=A.min (0.9)    $A.max (-0.9) x

    myfun1::Exp (Int,Int,Int)->Exp (Float,Float,Float)
    myfun1 xs=if (A.not (mypanduan xs) ) 
      then lift $ (0.0::Exp Float,0.0::Exp Float ,0.0::Exp Float)
      else 
        lift $ ( mybound cha1  +f1, mybound cha2+f2,mybound cha3+f3)
      where 
        (cha1,cha2,cha3)=myget_cha xs 
        (f1,f2,f3)= unlift $ myi2f xs::(Exp Float,Exp Float,Exp Float)

    final_fun::Exp (Int,Int,Int)->Exp (Float,Float,Float)
    final_fun xs=if   (A.not (mypanduan xs) ) A.||  A.abs cha1 A.>=0.5 A.|| A.abs cha2 A.>= 0.5 A.|| A.abs cha3 A.>= 0.5 
      then lift $ (0.0::Exp Float,0.0::Exp Float,0.0::Exp Float)  
      else if  A.abs (mytool xs1 xs2 xs3 +(dx1*cha1+dx2*cha2+dx3*cha3)*0.5) A.< 0.04/3.0 then 
        lift $ (0.0::Exp Float,0.0::Exp Float,0.0::Exp Float)
      else if (dd1+dd4)**2/(dd1*dd4-dd2*dd2) A.>12.1 then 
        lift $ (0.0::Exp Float,0.0::Exp Float,0.0::Exp Float)
      else lift $ (f1+cha1,f2+cha2,f3+cha3)
      where 
        (cha1,cha2,cha3)=myget_cha xs 
        (f1,f2,f3)= unlift $ myi2f xs::(Exp Float,Exp Float,Exp Float)
        (xs1,xs2,xs3)=unlift xs ::(Exp Int,Exp Int,Exp Int)
        (dx1,dx2,dx3)=myget_d  xs
        (dd1,dd2,dd3,dd4,dd5,dd6)=unlift$ myget_dd xs:: (Exp Float,Exp Float,Exp Float,Exp Float,Exp Float,Exp Float)

-- 
-- @output:: x y image_index m theta

compute_key_point_directions:: Acc (Array DIM3 Float)->Acc( Vector (Float,Float,Float))->Acc (Vector (Float,Float,Float,Float,Float))
compute_key_point_directions gauss_octave feas =final_temp
  where 
    k=2.0**(1/3.0)::Exp Float
    A.I1 n=A.shape feas
    A.I3 h w len=  A.shape gauss_octave 

    feas1::Acc (Vector (Int ,Float ,Float,Float))
    feas1= A.generate (A.I1 n) myfun  
      where 
        myfun ::Exp DIM1 ->Exp (Int ,Float ,Float ,Float)
        myfun (A.I1 i)= A.T4  i x y image_index
          where 
            A.T3 x y image_index= feas A.! (A.I1 i)

    compute_len::Exp (Int,Float,Float,Float)->Exp Int
    compute_len  (A.T4 _ _ _ image_index) =m*m
      where 
        sig=1.5*1.6*k**image_index
        m= 2*A.ceiling (  3*sig) +1  :: Exp Int

    compute_base_info::Exp (Int,Float,Float ,Float)->Exp Int-> Exp (Int,Float ,Float)
    compute_base_info (A.T4 i x y image_index) j= A.T3 i  temp_mag  (my_orientation ii1 ii2)
      where
        image_index1= A.round image_index ::Exp Int
        x1=A.round x ::Exp Int
        y1=A.round y ::Exp Int
        sig=1.5*1.6*k**image_index
        m= 2*A.ceiling (  3*sig) +1  :: Exp Int
        ii1= x1+ (A.mod j m)-A.ceiling(3*sig)
        ii2= y1+ (A.quot j m)-A.ceiling(3*sig)
        f_i1=A.fromIntegral $(A.mod j m) -A.ceiling(3*sig)
        f_i2=A.fromIntegral $(A.div j m) -A.ceiling(3*sig)
        
        temp_mag=(my_magnitude ii1 ii2)*A.exp (- (f_i1**2 +f_i2**2 ) /(2*sig*sig) ) /(2*A.pi*sig*sig)

        mytool::Exp Int->Exp Int->Exp Float
        mytool i1 i2 = gauss_octave A.! (A.I3  (mylegalx i1)  ( mylegaly i2)  image_index1 )

        mylegalx::Exp Int ->Exp Int
        mylegalx xx = A.max 0 $A.min (h-1) xx

        mylegaly::Exp Int->Exp Int
        mylegaly yy= A.max 0 $ A.min (w-1) yy

        my_magnitude::Exp Int->Exp Int->Exp Float
        my_magnitude i1 i2=if i1 A.<1 A.|| i1 A.>= (h-1) A.|| i2 A.<1 A.|| i2 A.>= (w-1) 
          then 
            0.0
          else 
            (dx*dx+dy*dy)**0.5
          where 
            dx=mytool (i1+1) i2 - mytool (i1-1) i2 
            dy=mytool i1 (i2+1) - mytool i1 (i2-1)

        my_orientation::Exp Int->Exp Int->Exp Float
        my_orientation i1 i2 =ori A.<0 ? (2*A.pi+ori,ori)
          where 
            ori= A.atan2 dy dx
            dx=mytool (i1+1) i2 - mytool (i1-1) i2 
            dy=mytool i1 (i2+1) - mytool i1 (i2-1)

    
    temp_baseinfo=A.expand compute_len compute_base_info feas1

    histogram_fun:: Acc (Vector (Int,Float,Float) )->Acc (Array DIM2 (Int,Float,Float))
    histogram_fun base_info = permute myfun1 chushi myfun2 base_info 
      where 
        chushi=A.generate (A.I2 n 36 ) (\(A.I2 i1 i2) -> lift (i1,0.0::Exp Float, 2.0*A.pi*(A.fromIntegral i2+0.5)/36.0  ::Exp Float)  )  
        
        myfun1::Exp (Int ,Float,Float)->Exp (Int,Float,Float)->Exp (Int,Float,Float)
        myfun1 (A.T3 _ mag1 _ ) (A.T3 i2 mag2 angle2) =A.T3 i2 (mag1+mag2) angle2
          
        myfun2:: Exp DIM1-> Exp (Maybe DIM2)
        myfun2  ix =A.Just_ (A.I2  i bini )
          where
            A.T3 i _ theta =base_info A.! ix

            bini= A.floor$ ( (theta) / A.pi)*0.5*36::Exp Int


    temp_his= histogram_fun temp_baseinfo
    smooth_his::Acc (Array DIM2 (Int,Float,Float))
    smooth_his = A.generate (A.I2 n 36) myfun
      where 
        myfun::Exp DIM2->Exp (Int, Float,Float)
        myfun (A.I2 i1 i2)  = lift $ ( i, 0.25*pre+0.5*mag +0.25*next , theta    ) 
          where 
            (A.T3 i mag theta)=temp_his A.! (A.I2 i1 i2)
            (A.T3 _ pre _) = temp_his A.! (A.I2 i1 (A.mod  (i2-1) 36 ) )
            (A.T3 _ next _)= temp_his A.! (A.I2 i1 (A.mod  (i2+1) 36 ) )


    max_direction= A.maximum $ A.map (\(A.T3 _ mag _ )-> mag) smooth_his

    peaks_his::Acc (Array DIM2 (Int,Float,Float))
    peaks_his= A.generate  (A.I2 n 36) myfun
      where 
        myfun::Exp DIM2->Exp (Int,Float,Float)
        myfun (A.I2 i1 i2) = if mag A.>= pre A.&& mag A.>=next 
          then 
            lift (i, mag,  theta+( 0.5*(pre-next)/(pre+next-2*mag))*A.pi/18.0   ) 
          else 
            lift (i,0.0::Exp Float,0.0::Exp Float)
          where 
            (A.T3 i mag theta)=smooth_his A.! (A.I2 i1 i2)
            (A.T3 _ pre _) = smooth_his A.! (A.I2 i1 (A.mod  (i2-1) 36 ) )
            (A.T3 _ next _)= smooth_his A.! (A.I2 i1 (A.mod  (i2+1) 36 ) )

--26025
--
    final_temp1=A.afst $A.filter (\(A.T3  i mag _)-> mag A.> 0.8*(max_direction A.! (A.I1 i)) ) peaks_his


    final_temp::Acc (Vector (Float,Float,Float,Float,Float))
    final_temp =A.flatten $A.map myfun final_temp1
      where 
        myfun ::Exp (Int,Float ,Float) ->Exp (Float ,Float, Float,Float ,Float)
        myfun (A.T3 i mag theta)= A.T5 x y image_index mag theta
          where 
            (A.T3 x y image_index)=feas A.! (A.I1 i)





compute_descriptors2::Acc (Array DIM3 Float)->Acc(Vector (Float ,Float ,Float,Float ,Float))->Acc( Array DIM2 Float)
compute_descriptors2 gauss_octave fea_dirs=normalize_des $A.map (A.min 0.2)  $ normalize_des$get_final_desc $ histogram_fun $ A.expand compute_len compute_base_info fea_dirs1 --  A.fill (A.I2 1 1) 0
  where 
    k=2.0**(1/3.0)::Exp Float
    n::Exp Int
    n=  A.min myn myn
      where 
        A.I1 myn=A.shape fea_dirs
    A.I3 h w len=  A.shape gauss_octave 
    step_angle =A.pi*2/8.0 ::Exp Float

    fea_dirs1::Acc (Vector (Int ,Float ,Float ,Float ,Float))
    fea_dirs1= A.generate (A.I1 n) myfun  
      where 
        myfun ::Exp DIM1 ->Exp (Int ,Float ,Float ,Float ,Float)
        myfun (A.I1 i)= A.T5  i x y image_index  theta
          where 
            A.T5 x y image_index _ theta= fea_dirs A.! (A.I1 i)

    compute_len::Exp (Int ,Float ,Float ,Float ,Float)->Exp Int
    compute_len  (A.T5 _ _ _ image_index  _) =(2*radius+1)*(2*radius+1)
      where 
        sig=1.5*1.6*k**image_index
        m= A.floor (  3*sig)  :: Exp Int
        radius=A.round(0.5*5*(A.fromIntegral m)::Exp Float ) ::Exp Int

    compute_base_info::Exp (Int ,Float ,Float ,Float ,Float)->Exp Int-> Exp (Int,Float,Float,Float,Float,Float,Float,Float,Float )
    compute_base_info (A.T5 i x y image_index angle) j= A.T9 i  (f_i4/f_m) (f_i5/f_m) temp_mag temp_angle f_i1 f_i2 f_i3 0
      where
        image_index1= A.round image_index ::Exp Int
        x1=A.round x ::Exp Int
        y1=A.round y ::Exp Int
        sig=1.5*1.6*k**image_index
        m= A.floor (  3*sig)  :: Exp Int
        radius=A.round(0.5*5*(A.fromIntegral m)::Exp Float ) ::Exp Int
        f_m= A.fromIntegral m::Exp Float

        ii2= A.mod j (2*radius+1)
        ii1= A.div j (2*radius+1)

        f_i1 =( A.min 3 $ A.max 0 $ A.fromIntegral$ (A.floor (f_i4/f_m+1.5)::Exp Int)  ) ::Exp Float
        
        f_i2 =(A.min 3 $ A.max 0 $A.fromIntegral$ (A.floor (f_i5/f_m+1.5)::Exp Int)  ) ::Exp Float
        f_i3 = (  A.min  7  $ A.fromIntegral$ (A.floor (temp_angle /step_angle)::Exp Int) )::Exp Float

        f_i4=A.fromIntegral ii1- (A.fromIntegral radius)::Exp Float
        f_i5=A.fromIntegral ii2- (A.fromIntegral radius)::Exp Float

        (temp_mag,temp_angle)=(A.exp (-(f_i4*f_i4+f_i5*f_i5)/(2*4*f_m*f_m))*my_magnitude1 f_i4 f_i5,((\x->x A.<0 ? (2*A.pi+x,x) ) $my_orientation1 f_i4 f_i5 - angle) )

        mytool::Exp Int->Exp Int->Exp Float
        mytool i1 i2 = gauss_octave A.! (A.I3  (mylegalx i1)  ( mylegaly i2)  image_index1 )

        mylegalx::Exp Int ->Exp Int
        mylegalx xx = A.max 0 $A.min (h-1) xx

        mylegaly::Exp Int->Exp Int
        mylegaly yy= A.max 0 $ A.min (w-1) yy

        my_magnitude::Exp Int->Exp Int->Exp Float
        my_magnitude i1 i2=if i1 A.<1 A.|| i1 A.>= (h-1) A.|| i2 A.<1 A.|| i2 A.>= (w-1) 
          then 
            0.0
          else 
            (dx*dx+dy*dy)**0.5
          where 
            dx=mytool (i1+1) i2 - mytool (i1-1) i2 
            dy=mytool i1 (i2+1) - mytool i1 (i2-1)

        my_orientation::Exp Int->Exp Int->Exp Float
        my_orientation i1 i2 =ori A.<0 ? (2*A.pi+ori,ori)
          where 
            ori= A.atan2 dy dx
            dx=mytool (i1+1) i2 - mytool (i1-1) i2 
            dy=mytool i1 (i2+1) - mytool i1 (i2-1)

-- 以上坐标是绝对值

        rotation_operator::Exp Float ->Exp Float->Exp (Float,Float)
        rotation_operator f1 f2= lift (f1*A.cos angle -f2*A.sin angle, f1*A.sin angle +f2*A.cos angle  )

        my_magnitude1::Exp Float->Exp Float->Exp Float
        my_magnitude1 f1 f2= my_magnitude (A.round n_f12_1 +x1) (A.round n_f12_2 +y1)
          where 
            (n_f12_1,n_f12_2)=unlift  $rotation_operator f1 f2

        my_orientation1::Exp Float->Exp Float->Exp Float
        my_orientation1 f1 f2 = my_orientation (A.round n_f12_1 +x1) (A.round n_f12_2 +y1)
          where 
            (n_f12_1,n_f12_2)=unlift  $rotation_operator f1 f2


-- (i j k) (i+1 j k) (i,j+1,k ) (i+1 ,j+1,k) (i j k+1) (i+1 j k+1) (i,j+1,k+1 ) (i+1 ,j+1,k+1)

    histogram_fun:: Acc (Vector (Int,Float,Float,Float,Float,Float,Float,Float,Float) )->Acc (Array DIM2 (Int,Float,Float,Float,Float,Float,Float,Float,Float))
    histogram_fun base_info = permute myfun1 chushi myfun2 base_info 
      where 
        i0=0::Exp Float
        chushi=A.generate (A.I2 n (4*4*8) ) (\(A.I2 i1 i2) ->lift (i1,i0,i0,i0,i0,i0,i0,i0,i0))

        myfun1::Exp (Int,Float,Float,Float,Float,Float,Float,Float,Float)->Exp (Int,Float,Float,Float,Float,Float,Float,Float,Float)->Exp (Int,Float,Float,Float,Float,Float,Float,Float,Float)
        myfun1 (A.T9 i fx fy mag theta fi1 fi2 fi3 _ ) (A.T9 _ ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8) =A.T9 i (ma1+mag1) (ma2+mag2) (ma3+mag3) (ma4+mag4) (ma5+mag5) (ma6+mag6) (ma7+mag7) (ma8+mag8)
          where
            myfun3 :: Exp Float ->Exp Float ->Exp Float->Exp Float
            myfun3 f1 f2 f3= dx*dy*dz*mag
              where 
                f11=   f1-1.5
                f22=   f2-1.5
                f33=   (f3+0.5)*step_angle
                dx= 1- ( A.min  1  (A.abs (fx-f11) ))
                dy= 1- ( A.min  1  (A.abs (fy-f22) ))
                dz= step_angle - (A.min step_angle  $ (\x->x A.>A.pi ? ((2*A.pi-x),x)  ) $ A.abs (theta  - f33)  )

            mag1=myfun3 (fi1) (fi2) (fi3)
            mag2=myfun3 (fi1+1) (fi2) (fi3)
            mag3=myfun3 (fi1) (fi2+1) (fi3)
            mag4=myfun3 (fi1+1) (fi2+1) (fi3)
            mag5=myfun3 (fi1) (fi2) (fi3+1)
            mag6=myfun3 (fi1+1) (fi2) (fi3+1)
            mag7=myfun3 (fi1) (fi2+1) (fi3+1)
            mag8=myfun3 (fi1+1) (fi2+1) (fi3+1)

        myfun2:: Exp DIM1-> Exp (Maybe DIM2)
        myfun2  ix =A.Just_ (A.I2  i  ((A.round fi1)*4*8+(A.round fi2)*8+(A.round fi3) )   )
          where
            A.T9 i _ _ _ _ fi1 fi2 fi3 _  =base_info A.! ix
           

    get_final_desc::Acc (Array DIM2 (Int,Float,Float,Float,Float,Float,Float,Float,Float))->Acc (Array DIM2 Float)
    get_final_desc his= A.generate (A.I2 n 128) myfun 
      where 
        myfun::Exp DIM2->Exp Float
        myfun (A.I2 i1 i2)= mag1+mag2+mag3+mag4+mag5+mag6+mag7+mag8
          where 
            i21= A.quot i2 32
            i22=A.mod (A.quot i2 8) 4
            i23=A.mod i2 8
            mytool1::Exp Int->Exp Int->Exp Int->Exp  Float
            mytool1 j1 j2 j3=(i21 A.== j1 A.&& i22 A.== j2 A.&& i23 A.== j3) ? (1::Exp Float,0 ::Exp Float)

            mytool::Exp Int->Exp Int->Exp Int->Exp Float
            mytool j1 j2 j3 =(mytool1 j1 j2 j3) * ma1+
              (mytool1 (j1+1) j2 j3) *ma2 +
              (mytool1 j1 (j2+1) j3) *ma3+
              (mytool1 (j1+1) (j2+1) j3) *ma4+
              (mytool1 j1 j2 (j3+1) ) *ma5+
              (mytool1 (j1+1) j2 (j3+1)  ) *ma6+
              (mytool1 j1 (j2+1) (j3+1) ) *ma7+
              (mytool1 (j1+1) (j2+1) (j3+1) ) *ma8
              where 
                (A.T9 _ ma1 ma2 ma3 ma4 ma5 ma6 ma7 ma8 )= his A.! (A.I2 i1 (j1*4*8+j2*8+j3 )) 
            mylegal1::Exp Int ->Exp Int
            mylegal1 x=A.min 3 $A.max 0 x 

            mylegal2::Exp Int ->Exp Int
            mylegal2 x=A.min 7 $A.max 0 x 

            mytool2::Exp Int->Exp Int->Exp Int->Exp Float
            mytool2 j1 j2 j3 =mytool (mylegal1 j1) (mylegal1 j2) (mylegal2 j3)

            mag1=mytool2 i21 i22 i23
            mag2=mytool2 (i21-1) i22 i23
            mag3=mytool2 i21 (i22-1) i23
            mag4=mytool2 (i21-1) (i22-1) i23
            mag5=mytool2 i21 i22 (i23-1)
            mag6=mytool2 (i21-1) i22 (i23-1)
            mag7=mytool2 i21 (i22-1) (i23-1)
            mag8=mytool2 (i21-1) (i22-1) (i23-1)


    normalize_des:: Acc (Array DIM2 Float)->Acc (Array DIM2 Float)
    normalize_des des=A.generate (A.I2 n 128) mult_fun
      where 
        norm_v=A.map (\x->x**0.5) ( A.fold (+) 0 $ (A.zipWith (*) des des)) 
        mult_fun::Exp DIM2 ->Exp Float
        mult_fun (A.I2 i1 i2)= (des A.! (A.I2 i1 i2) )/(norm_v A.! (A.I1 i1) ) 
  


compute_descriptors::Acc (Array DIM3 Float)->Acc(Vector (Float ,Float ,Float,Float ,Float))->Acc( Array DIM2 Float)
compute_descriptors gauss_octave fea_dirs= A.afst loop  
  where
    A.I1 n=A.shape fea_dirs
    A.I3 h w len=  A.shape gauss_octave

    initial_ =(use $ A.fromList (Z:.131 :. 0) [] )::Acc (Array DIM2 Float)
    loop= A.awhile 
      (\(A.T2 _ s)-> A.map (A.< n) s ) 
      (\(A.T2 f s )->
        let 
          ss=the s
          new1=A.flatten$ compute_descriptor gauss_octave (fea_dirs A.! (A.I1 ss))
          (x,y ,image_index,mag,angle)= unlift$ fea_dirs A.! (A.I1 ss)::(Exp Float,Exp Float,Exp Float,Exp Float,Exp Float)
          new2= (A.fill (A.I1 1) (x/(A.fromIntegral h))) A.++ (A.fill (A.I1 1) (y/(A.fromIntegral w))) A.++ (A.fill (A.I1 1) (angle))
          new= A.reshape (A.I2 131 1) $ new2 A.++ new1
        in 
          A.T2  ( f A.++ new )  (unit (ss+1))  ) 
      (A.T2  initial_ (A.unit 0))


compute_descriptor:: Acc (Array DIM3 Float)->Exp (Float,Float,Float,Float,Float)->Acc (Array DIM3 Float)
compute_descriptor gauss_octave fea_dir=  final_re
  where 
    (x,y,image_index,mag,angle)=unlift fea_dir::(Exp Float,Exp Float,Exp Float,Exp Float,Exp Float)
    image_index1= A.round image_index::Exp Int
    x1=A.round x ::Exp Int
    y1=A.round y ::Exp Int
    k=2.0**(1/3.0)::Exp Float
    sig=1.6*k**image_index
    m= A.ceiling (  3*sig)  :: Exp Int

    radius=A.round(0.5*5*2.0**0.5*(A.fromIntegral m)::Exp Float ) ::Exp Int

    A.I3 h w len=  A.shape gauss_octave
-- 以下坐标是绝对值

    mytool::Exp Int->Exp Int->Exp Float
    mytool i1 i2 = gauss_octave A.! (A.I3  (mylegalx i1)  ( mylegaly i2)  image_index1 )

    mylegalx::Exp Int ->Exp Int
    mylegalx xx = A.max 0 $A.min (h-1) xx

    mylegaly::Exp Int->Exp Int
    mylegaly yy= A.max 0 $ A.min (w-1) yy
 
    my_magnitude::Exp Int->Exp Int->Exp Float
    my_magnitude i1 i2=if i1 A.<1 A.|| i1 A.>= (h-1) A.|| i2 A.<1 A.|| i2 A.>= (w-1) 
      then 
        0.0
      else 
        (dx*dx+dy*dy)**0.5
      where 
        dx=mytool (i1+1) i2 - mytool (i1-1) i2 
        dy=mytool i1 (i2+1) - mytool i1 (i2-1)

    my_orientation::Exp Int->Exp Int->Exp Float
    my_orientation i1 i2 = if ori A.<0 
      then 
        2*A.pi +ori
      else 
        ori
      where 
        ori=A.atan2 dy dx
        dx=mytool (i1+1) i2 - mytool (i1-1) i2 
        dy=mytool i1 (i2+1) - mytool i1 (i2-1)

-- 以上坐标是绝对值

    rotation_operator::Exp Float ->Exp Float->Exp (Float,Float)
    rotation_operator f1 f2= lift (f1*A.cos angle -f2*A.sin angle, f1*A.sin angle +f2*A.cos angle  )

    my_magnitude1::Exp Float->Exp Float->Exp Float
    my_magnitude1 f1 f2= my_magnitude (A.round n_f12_1 +x1) (A.round n_f12_2 +y1)
      where 
        (n_f12_1,n_f12_2)=unlift  $rotation_operator f1 f2


    my_orientation1::Exp Float->Exp Float->Exp Float
    my_orientation1 f1 f2 = my_orientation (A.round n_f12_1 +x1) (A.round n_f12_2 +y1)
      where 
        (n_f12_1,n_f12_2)=unlift  $rotation_operator f1 f2


    myfun::Exp DIM5->Exp Float
    myfun (A.I5  i1 i2 i3 i4 i5)=dx*dy*dz*temp_mag
      where 
        f_m= A.fromIntegral m::Exp Float
        step_angle =A.pi*2/8.0 ::Exp Float

        f_i1=(A.fromIntegral i1-1.5)*f_m::Exp Float
        f_i2=(A.fromIntegral i2 -1.5)*f_m ::Exp Float 
        f_i3=(A.fromIntegral i3 +0.5)*step_angle ::Exp Float

        f_i4=A.fromIntegral i4- (A.fromIntegral radius)::Exp Float
        f_i5=A.fromIntegral i5- (A.fromIntegral radius)::Exp Float

        dx=f_m  - (A.min f_m $ A.abs (f_i4-f_i1)) ::Exp Float
        dy=f_m - (A.min (f_m)$ A.abs (f_i5-f_i2))::Exp Float
        dz= step_angle - (A.min step_angle  $ A.abs (temp_angle  - f_i3))
        temp_angle::Exp Float
        temp_angle= if temp_a1 A.<0 
          then 
            temp_a1 + 2*A.pi
          else  
            temp_a1
          where
            temp_a1= my_orientation1 f_i4 f_i5 - angle 

        temp_mag=A.exp (-(f_i4*f_i4+f_i5*f_i5)/(2*4*f_m*f_m))*my_magnitude1 f_i4 f_i5

    chushi::Acc (Array DIM5 Float )
    chushi = A.generate (A.I5 4 4 8 (2*radius+1) (2*radius+1) ) myfun

    final_re::Acc (Array DIM3 Float)
    final_re=normalize_des $A.map  (A.min 0.2 ) $  normalize_des  temp_re
      where 
        temp_re=A.fold (+) 0.0  $ A.reshape (A.I4 4 4 8 ((2*radius+1)*(2*radius+1))) chushi 
        normalize_des:: Acc (Array DIM3 Float)->Acc (Array DIM3 Float)
        normalize_des des=A.map mult_fun des
          where 
            norm_v=A.the $ A.fold1All (+) (A.zipWith (*) des des) :: Exp Float
            mult_fun::Exp Float ->Exp Float
            mult_fun xx= xx/(norm_v**0.5) 


pixelrgb2gray:: Acc (Array DIM2 PixelRGB8)->Acc (Array DIM2 Float)
pixelrgb2gray image= A.map pix2double image
  where 
    pix2double::Exp PixelRGB8->Exp Float
    pix2double (PixelRGB8_ a1 a2 a3)= (A.fromIntegral a1* 0.299+ A.fromIntegral a2*0.587+A.fromIntegral a3*0.114)/255.0


gray2grayimage:: Acc (Array DIM2 Float)->Acc (Array DIM2  PixelF)
gray2grayimage  image=A.map (\x-> x/255.0) image  



normalization_key_points::Acc (Array DIM3 Float)-> Acc (Vector (Float,Float,Float,Float,Float)) ->Acc (Vector (Float,Float,Float,Float,Float))
normalization_key_points gauss_octave key_dirs =A.map ( \(A.T5 x y image_index mag theta) ->A.T5 (x/fh) (y/fw) image_index mag theta    ) key_dirs
  where 
    A.I3 h w len=  A.shape gauss_octave
    fh=A.fromIntegral h 
    fw=A.fromIntegral w




match_descr_points::Acc( Array DIM2 Float)->Acc( Array DIM2 Float)->Acc( Vector (Int,Int,Float))
match_descr_points des1 des2= A.afst $ A.filter (\(A.T3 _ _ vv)->vv A.== 1.0)  $ histogram_fun dotsijv
  where 
    A.I2 n1 m1=A.shape des1 
    A.I2 n2 m2 =A.shape des2
  
    my_dot::Exp DIM2->Exp Float
    my_dot (A.I2 x y) = Prelude.foldl (A.+) (0::Exp Float)  $ Prelude.map (\i ->(des1 A.! (A.I2 x i ) )*(des2 A.! (A.I2 y i)) ) $ (Prelude.map (A.constant) ([0..127]::[Int]) )  

    dots::Acc(Array DIM2 Float)
    dots = A.generate (A.I2 n1 n2)   (A.abs.my_dot)

    dotsijv= A.generate (A.I2 n1 n2) (\(A.I2 i j)-> lift (i,j,0::Exp Float ) ) ::Acc (Array DIM2 (Int,Int,Float)) 

    max_dots= A.maximum dots

    histogram_fun:: Acc (Array DIM2 (Int,Int,Float) )->Acc( Vector (Int,Int,Float))
    histogram_fun base_info= permute myfun1 chushi myfun2 base_info
      where 
        chushi= A.generate (A.I1 n1) (\(A.I1 i )->A.T3 i 0 0 )

        myfun1::Exp(Int,Int,Float)->Exp (Int,Int,Float)->Exp (Int,Int,Float)
        myfun1 (A.T3  i1 j1 v1) (A.T3 i2 j2 v2)= (dots A.! (A.I2 i1 j1)) A.> 0.6*(max_dots A.! (A.I1 i1)) A.? ( A.T3 i1 j1 (v2+1) , A.T3 i2 j2 v2 ) 
        
        myfun2:: Exp DIM2-> Exp (Maybe DIM1)
        myfun2  ix =A.Just_ (A.I1 i1 )
          where
            A.T3 i1 _ _ =base_info A.! ix






myget_curtime:: IO UTCTime
myget_curtime = do 
    a<- getCurrentTime
    return $ addUTCTime (60*60*8) a



myget_curdaytime::IO DiffTime
myget_curdaytime= return.utctDayTime=<<myget_curtime


difftimtofloat:: DiffTime->Float
difftimtofloat a=realToFrac a 

inttofloat:: Int->Float
inttofloat a=realToFrac a    


compute_keys_and_decrs:: String -> IO (Vector (Float,Float,Float,Float,Float),Array DIM2 Float)
compute_keys_and_decrs image_path= do 
  let 
    my_compute_accurate_extr:: Acc (Array DIM3 Float)->Acc (Vector (Float,Float,Float))
    my_compute_accurate_extr dog = compute_accurate_extrum_pos dog  $ compute_extremum_points dog 
    my_compute_keys::Acc(Array DIM3 Float)->Acc (Array DIM3 Float)->Acc (Vector (Float,Float,Float,Float,Float))
    my_compute_keys dog gauss = compute_key_point_directions gauss  $ my_compute_accurate_extr dog

    myget_one_re::Array DIM3 Float->Array DIM3 Float->(Vector (Float,Float,Float,Float,Float),Array DIM2 Float )
    myget_one_re dog gauss =(f_one1, s_one)
      where 
        f_one=CPU.run $ my_compute_keys (use dog) (use gauss)
        f_one1=CPU.run $ normalization_key_points (use dog) (use f_one)
        s_one=CPU.run $ compute_descriptors2 (use gauss) (use f_one)

    myfun::(Vector (Float,Float,Float,Float,Float),Array DIM2 Float )->(Vector (Float,Float,Float,Float,Float),Array DIM2 Float )-> (Vector (Float,Float,Float,Float,Float),Array DIM2 Float ) 
    myfun (f1,s1) (f2,s2)=  (rf , rs    ) 
      where 
        rf= CPU.run $(use f1) A.++ (use f2)
        -- rs= CPU.run $A.transpose $ (A.transpose $use s1) A.++ (A.transpose $use s2) 
        rs= CPU.run $ A.concatOn ACL._2 (use s1) (use s2 )
  h1<- openBinaryFile image_path ReadMode
  bstr<- B.hGetContents h1
  jpg_data<-return $decodeJpeg bstr

  rgb_data<- return $ convertRGB8  $ (\(Right x)->x) jpg_data

  n_octs<- return $Prelude.ceiling $logBase 2  (Prelude.fromIntegral $Prelude.min  (imageHeight rgb_data) (imageWidth rgb_data)) -2.0

  putStrLn $ show n_octs
  -- scale_image <- return $ CPU.run $ 

  base_gray_image<- return $CPU.run  $compute_base_gray_image$  pixelrgb2gray$ scale_image (use $ arrayOfImage rgb_data) (A.constant $ 2*imageWidth rgb_data) (A.constant $ 2*imageHeight rgb_data) 
  
  all_octaves <- return  $ generate_gauss_pyramid   base_gray_image n_octs
  putStrLn "begin dog pyramid"

  all_dog_octaves <- return $generate_dog_pyramid all_octaves

  re<-return $  Prelude.foldl1 myfun$ Prelude.map (\i->myget_one_re  (all_dog_octaves Prelude.!! i) (all_octaves Prelude.!! i ) ) [0..(n_octs-2)]

  -- re<-return $ myget_one_re  (all_dog_octaves Prelude.!! 5) (all_octaves Prelude.!! 5 )



  hClose h1
  return re



main :: IO ()
main = do
  putStrLn "Hello, Haskell!"
  putStrLn .show .difftimtofloat=<<myget_curdaytime  

  (keys,descs)<-  compute_keys_and_decrs "002_DSC.JPG"
  putStrLn $ show $CPU.run $unit $ A.shape (use descs) 
  putStrLn .show .difftimtofloat=<<myget_curdaytime  

  putStrLn $ show $CPU.run $unit $ A.shape (use keys) 


  -- h1<- openBinaryFile "001_DSC.JPG" ReadMode
  -- bstr<- B.hGetContents h1
  -- jpg_data<-return $decodeJpeg bstr

  -- rgb_data<- return $ convertRGB8  $ (\(Right x)->x) jpg_data


  -- putStrLn .show .difftimtofloat=<<myget_curdaytime  

  -- scale_image <- return $ CPU.run $ scale_image (use $ arrayOfImage rgb_data) (A.constant $ 2*imageWidth rgb_data) (A.constant $ 2*imageHeight rgb_data) 
    
  -- putStrLn .show .difftimtofloat=<<myget_curdaytime  

  -- gray_image <- return$CPU.run$ pixelrgb2gray $ ( use$ scale_image)  
  -- putStrLn .show .difftimtofloat =<< myget_curdaytime  

  -- base_gray_image <- return $ CPU.run $ compute_base_gray_image$ use gray_image 
  -- putStrLn "begin gauss pyramid"



  -- all_octaves <- return  $ generate_gauss_pyramid   base_gray_image 10
  -- putStrLn "begin dog pyramid"
  -- all_dog_octaves <- return $generate_dog_pyramid all_octaves

  -- putStrLn "begin extremum points"

  -- extremum_points <- return $ CPU.run $ compute_extremum_points $ use $ all_dog_octaves Prelude.!! 9
  -- putStrLn "begin accurate extrum pos"
  
  -- accurate_extrum_pos <- return $ CPU.run$ compute_accurate_extrum_pos (use $all_dog_octaves Prelude.!! 9) (use extremum_points)
   
  
  -- putStrLn $ show $CPU.run $unit $ A.shape (use accurate_extrum_pos)

  -- key_dirs<-return $ CPU.run$ compute_key_point_directions (use $ all_octaves Prelude.!! 9 ) (use accurate_extrum_pos)

  -- descrs<- return $ CPU.run $ compute_descriptors2 (use $ all_octaves Prelude.!! 9 ) (use key_dirs)
    
  -- putStrLn $ show $ CPU.run $unit $ A.shape (use descrs)
  
  -- putStrLn .show .difftimtofloat=<<myget_curdaytime  


  MyLib.someFunc

